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Abstract 

Background: Pleistocene climatic oscillations have played a major role in structuring present-day biodiversity. The 
southern Mediterranean peninsulas have long been recognized as major glacial refugia, from where Northern 
Europe was post-glacially colonized. However, recent studies have unravelled numerous additional refugia also in 
northern regions. We investigated the phylogeographic pattern of the widespread Western Palaearctic lizard 
Podarcis muralis, using a range-wide multilocus approach, to evaluate whether it is concordant with a recent 
expansion from southern glacial refugia or alternatively from a combination of Mediterranean and northern refugia. 

Results: We analyzed DNA sequences of two mitochondrial {cytb and nd4) and three nuclear {acm4, mclr, and pdc) 
gene fragments in individuals from 52 localities across the species range, using phylogenetic and phylogeographic 
methods. The complex phylogeographic pattern observed, with 23 reciprocally monophyletic alio- parapatric 
lineages having a Pleistocene divergence, suggests a scenario of long-term isolation in multiple ice-age refugia 
across the species distribution range. Multiple lineages were identified within the three Mediterranean peninsulas - 
Iberia, Italy and the Balkans - where the highest genetic diversity was observed. Such an unprecedented 
phylogeographic pattern - here called "refugia within all refugia" - compasses the classical scenario of multiple 
southern refugia. However, unlike the southern refugia model, various distinct lineages were also found in northern 
regions, suggesting that additional refugia in France, Northern Italy, Eastern Alps and Central Balkans allowed the 
long-term persistence of this species throughout Pleistocene glaciations. 

Conclusions: The phylogeography of Podarcis muralis provides a paradigm of temperate species survival in 
Mediterranean and extra-Mediterranean glacial refugia. Such refugia acted as independent biogeographic 
compartments for the long-term persistence of this species, for the differentiation of its genetic lineages, and for 
the short-distance post-glacial re-colonization of neighbouring areas. This finding echoes previous findings from 
recent phylogeographic studies on species from temperate ecoregions, thus suggesting the need for a reappraisal 
of the role of northern refugia for glacial persistence and post-glacial assembly of Holarctic biota. 
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Background 

Pleistocene climatic oscillations have played a major role 
in structuring current biodiversity patterns and shaping 
the distribution of species and their genetic diversity 
both at the regional and global scale [1,2]. In recent de- 
cades, phylogeographic studies have been used to assess 
the genetic consequences of Pleistocene ice ages on 
various organisms, highlighting the dynamic nature of 
species ranges and the role of microevolutionary pro- 
cesses in determining the extent and structure of intra- 
specific diversity [1,3,4]. 

In Europe, the genetic structure of temperate organisms 
has been mainly explained by cycles of contraction toward 
glacial refugia located in southern peninsulas and post- 
glacial recolonization of northern regions, tracking the ex- 
pansion of suitable habitats [5,6] Accordingly, southern 
peninsular refugia are thought to have played a major role 
in the long-term maintenance of genetic diversity and dif- 
ferentiation, and relatively few patterns and routes of 
colonization have been described as paradigms for the 
post-glacial arrival of species in Northern Europe [6-8]. 

However, in recent years, multiple evidence from pa- 
leoclimatic, palaeontological and phylogeographic stud- 
ies have identified the existence of glacial refugia in 
northern regions of Europe, or "northern refugia", and 
have highlighted the prominent contribution of these re- 
fugia to the present-day genetic diversity of many organ- 
isms [9-12]. Mounting examples from plants, insects, 
mammals, amphibians, and reptiles show that in many 
temperate species the southern refugia model is insufficient 
to explain the observed genetic patterns, which would be 
better explained by a scenario of ice-age survival in a 
combination of Mediterranean and extra-Mediterranean 
refugia [13-17]. This changing view with respect to geo- 
graphical location of refugia has also had direct implica- 
tions for the inferred scenarios of recolonization for many 
widespread species from the Western Palaearctic. This is 
true even in emblematic species, such as the brown bear 
and the common beech, the phylogeographic patterns 
of which were previously upheld as paradigms of post- 
glacial recolonization routes (exclusively) from southern 
refugia [6,7]. In these species, the existence of additional 
extra-Mediterranean refugia has recently become evident 
through the analysis of fossil DNA and palaeobotanical 
data, respectively. These evidence suggests that range 
shifts and expansions in these species could have been 
much smaller than previously thought [18-21]. 

Detailed phylogeographic data from widespread Western 
Palaearctic species are particularly valuable for evaluat- 
ing the plausibility of a scenario of ice-age survival in 
Mediterranean and extra-Mediterranean refugia. The 
common Wall lizard, Podarcis muralis (Laurenti, 1768) 
is a Western Palaearctic species which provides a suit- 
able case-study for this purpose. Podarcis muralis is a 



locally abundant lacertid lizard found in a variety of en- 
vironments across a wide altitudinal range (from sea 
level to over 2000 m of altitude). It exhibits considerable 
variation in colour patterns, biometry and pholidosis, 
which has led to the description of several subspecies 
[22,23]. The distribution of P. muralis is unusual rela- 
tive to that of its congeners and, together with other 
characteristics, makes this species a useful model for 
evaluating the relative contribution of southern versus 
extra-Mediterranean refugia in shaping the current dis- 
tribution of species and their genetic diversity. Firstly, it 
covers a wide range across the Western Palaearctic- 
from North-Eastern Anatolia and the Black Sea coast 
through all the Mediterranean peninsulas - Iberia, Italy, 
the Balkans- as well as across extra-Mediterranean re- 
gions of France, Germany, Slovenia, Austria and the 
Central Balkans (Figure 1). Secondly, it is abundant in 
both Mediterranean and continental climatic conditions 
and its distribution does not seem to be strongly limited 
by the presence of other lizards. Throughout the southern 
peninsulas P. muralis occurs in sympatry with more local- 
ized Podarcis species, such as members of the P. hispanica 
group in the Iberian Peninsula, P. sicula in the Italian 
Peninsula, and P. erhardii in the Balkans as well as with 
members of the genera Lacerta, Iberolacerta, and Timon 
across its whole range. Also, it is a highly successful inva- 
sive species in North- Western Europe where around 150 
non-native P. muralis populations have been identified 
[24]. Thus, Mediterranean and extra-Mediterranean areas 
seem generally suitable for the survival and colonization 
of this species. Moreover, although little is presently 
known about the evolutionary history of P. muralis, either 
a scenario of recent colonization of its northern range 
from a southern primary range or a long-term history of 
the species in some of these northern areas can be in- 
voked to explain the observed intraspecific variation. 
Harris and Arnold [25] hypothesized that P. muralis re- 
cently colonized most of its distribution from the Italian 
Peninsula, based on the observation that it exhibits the 
greatest morphological and allozyme diversity in Italy [26] . 
On the other hand, recent genetic assessment identified 
divergent genetic lineages - likely of Pleistocene origin - in 
Italy but also in France and the Balkans, disproving a 
colonization of these two latter areas from Italy [24,27,28] . 
Nevertheless, these studies have been limited to small 
parts of the range and analyzed only mitochondrial DNA, 
and particularly partial cytochrome b sequences, thus 
preventing a full appreciation of the origin of these line- 
ages and of the overall evolutionary history and biogeog- 
raphy of the species. 

In this study, we analysed DNA sequences data from 
86 individuals from 52 localities across the species range, 
including extensive sampling in the Mediterranean pen- 
insulas where native diversity is expected to be highest. 
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(See figure on previous page.) 

Figure 1 Sampling localities, phylogenetic relationships, and distribution of mitochondrial haplotypes and haploclades. Maximum 

Likelihood (ML) trees depicting the phylogenetic relationships among 75 combined mitochondrial sequences (cytb and nd4) of Podarcis murolis 

from 52 localities (A) and among 185 cytb sequences of Podorcis murolis obtained in this study (N=82) and from GenBank, (N=1 03) (see 

Additional file 1) (B). ML bootstrap values over 1000 replicates are reported in correspondence to the nodes; black stars indicate BA posterior 

probabilities = 1.00; Podorcis liolepis was used as an outgroup. In the combined cytb/nd4 ML tree specimen codes are reported (see Table 1 for 

details) and main mitochondrial clades numbered from 1 to 17; in the cytb ML tree tip nodes are collapsed and main mitochondrial clades are 

numbered from 1 to 23 according to the combined cytb/nd4 ML tree and named according to previous studies. Samples used in this study and 

geographic distribution of mitochondrial clades are reported in the map: big circles show the 53 localities sampled in this study and triangles 

show the geographic origin for the cytb sequences obtained from GenBank, coloured according to mitochondrial clades defined by the ML trees 

(black triangles indicate admixed populations for which the frequency of haplotypes belonging to each clades is shown by the pie diagrams; 

grey circles indicate samples for which only nuclear sequences were available). 
I J 



We used two partial mitochondrial genes (cytb and nd4), 
and three partial nuclear genes {acm4, mclr, and pdc), 
all of which have been successfully used for assessing 
variation within other species of lacertids (e.g. [29,30]). 
Furthermore, by comparing 103 published cytochrome b 
sequences from previous studies with our assessment, 
we placed these in a more robust phylogenetic frame- 
work, and further used them to confirm that all known 
diversity was included in this study. This comprehensive 
data set enables us to evaluate whether the patterns of 
genetic variation observed in this species reflect a recent 
expansion from a single or a few Pleistocene southern 
glacial refugia or, alternatively, from a combination of 
Mediterranean and extra-Mediterranean refugia. Under 
the first scenario we expect to identify endemic lineages in 
one or more southern refugia, to observe only southern- 
derived haplotypes in northern areas, and an overall pat- 
tern of southern richness vs northern purity' of genetic 
diversity [1,2,6]. Alternatively, if P. muralis also survived 
in northern glacial refugia, we predict an occurrence of 
lineages endemic to northern regions, with an ancient 
(Pleistocene) divergence from southern lineages, and an 
absence of a northern purity pattern of genetic diversity. 
The main goal of our range-wide multilocus approach is, 
therefore, to identify the biogeographic and evolutionary 
processes shaping the genetic diversity of P. muralis and 
to discuss, in the framework of European phylogeography, 
how this species came to have such an unusual biogeo- 
graphic pattern. 

Methods 

Data collection 

We sampled a total of 86 individuals from 52 localities 
spanning the range of Podarcis muralis with particular 
attention to the Iberian, Italian, and Balkan peninsulas. 
Three specimens of two Iberian Podarcis, P. liolepis 
(Boulenger, 1905) and P. bocagei (Seoane, 1884), and the 
Italian wall lizard P. sicula (Rafinesque-Schmaltz, 1810), 
were sampled and designated as outgroups. We used 
multiple outgroups because the phylogenetic position of 
P. muralis within the genus Podarcis is not well resolved, 
yet this species seems more closely related to Iberian 



Podarcis or P. sicula [31,32]. Additionally, partial cytb se- 
quences for 103 individuals were obtained from GenBank 
and included in phylogenetic analyses. Specimens included 
in molecular analyses, together with collection locality de- 
tails, are reported in Table 1 and Figure 1. Accession num- 
bers for sequences retrieved from GenBank are reported 
in Additional file 1. 

Total genomic DNA was extracted from alcohol- 
preserved tail muscle collected from live specimens follow- 
ing standard high-salt protocols [33]. Portions of five genes 
were amplified by polymerase chain reaction: two mito- 
chondrial genes, cytochrome b {cytb) and NADH dehydro- 
genase subunit 4 with flanking tRNAs (yid4), and three 
nuclear genes, Melanocortin receptor 1 (mclr), acetylcho- 
linergic receptor M4 (acm4), and Phosducin (pdc). The 
primers used for amplification and respective references are 
reported in Additional file 2. For pdc we used primers 
newly developed for this study. Amplifications were carried 
out in 25 \iL volumes, containing IX PCR buffer (50 mm 
Tris-HCl, 50 mm NaCl, pH 8.5); 3 mM MgCl 2 ; 0.6 mM 
each dNTP, 2U of GoTaq DNA polymerase (Promega), 
0.4 \iM each primer and approximately 50 ng of genomic 
DNA. Amplification conditions consisted of a preliminary 
denaturation step at 92°C for 3 minutes, followed by 16 
touchdown cycles with 30 seconds at 92°C, annealing 
temperature decreasing 0.5°C per cycle from 60°C to 52°C 
(30 seconds) and extension for 1 minute at 72°C. 20 more 
cycles similar to these but with annealing temperatures 
stable at 52°C followed. A final extension was carried out at 
72°C for 15 minutes. 

Data analysis 

Electropherograms were checked and consensus sequences 
were aligned in GENEIOUS 6.0 (www.geneious.com) using 
the Geneious Alignment algorithm. We discarded the hy- 
pothesis of pseudogenes occurring in our mitochondrial 
sequence dataset by confirming: (i) the absence of stop co- 
dons in protein-coding fragments, (ii) the consistency of 
the base composition pattern with a mitochondrial origin 
(less than 5% of guanines in the third position, see [34]), 
and (iii) the overall similarity with the reference mitochon- 
drial genome of P. muralis (Genbank accession number 
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Table 1 Geographical location and codes for the studied specimens of Podarcis muralis and the outgroup species 



P. liolepis 


Specimen 


Locality 


Locality 








Genbank accession numbers 




code 


code 


Coordinates 


cytb 


nd4 


mclr 


acm4 


pdc 


DB 16840 


1 


Vieille-Roche (France) 


47.50 N 


-2.38 E 


KF372191 


KF372350 


KF372113 


KF372273 


KF372036 


DB16841 


1 


Vieille-Roche (France) 


47.50 N 


-2.38 E 


KF372192 


KF372351 


KF372114 


KF372274 


KF372037 


DB4296 


2 


Palacios de Compludo (Spain) 


42.46 N 


-6.47 E 


KF372193 


KF372352 


KF3721 15 


KF372275 


KF372038 


DB15980 


3 


Braha de Sosas (Spain) 


43.00 N 


-6.31 E 


KF372194 


KF372353 


KF372116 


KF372276 


KF372039 


DB15977 


4 


Sehora de Carrasconte (Spain) 


42.94 N 


-6.23 E 


KF372195 


KF372354 


KF372117 


KF372277 


KF372040 


DB4281 


5 


La Omahuela (Spain) 


42.78 N 


-5.98 E 


KF372196 


KF372355 


KF372118 


- 


- 


DB4295 


6 


Leon (Spain) 


42.59 N 


-5.58 E 


KF372197 


KF372356 


KF372119 


KF372278 


KF372041 


DB15968 


7 


La Candamia (Spain) 


42.60 N 


-5.55 E 


KF372198 


KF372357 


KF372120 


KF372279 


KF372042 


DB15969 


8 


Valdehuesa (Spain) 


42.94 N 


-5.32 E 


KF372199 


KF372358 


KF372121 


KF372280 


KF372043 


DB15970 


8 


Valdehuesa (Spain) 


42.94 N 


-5.32 E 


KF372200 


KF372359 


KF372122 


KF372281 


KF372044 


DB8970 


9 


Tanes (Spain) 


43.21 N 


-5.40 E 


KF372201 


- 


KF372123 


KF372282 


KF372045 


DB8971 


9 


Tanes (Spain) 


43.21 N 


-5.40 E 


KF372202 


KF372360 


KF372124 


KF372283 


KF372046 


DB7092 


10 


Turieno (Spain) 


43.16 N 


-4.66 E 


- 


KF372361 


KF372125 


KF372284 


KF372047 


DB5880 


11 


Matienzo (Spain) 


43.33 N 


-3.59 E 


KF372203 


KF372362 


KF372126 


KF372285 


KF372048 


DB5882 


11 


Matienzo (Spain) 


43.33 N 


-3.59 E 


KF372204 


KF372363 


KF372127 


KF372286 


KF372049 


DB7065 


12 


Orihon (Spain) 


43.40 N 


-3.33 E 


KF372205 


KF372364 


KF372128 


KF372287 


KF372050 


DB8957 


13 


Guadarrama (Spain) 


40.71 N 


-4.14 E 


KF372206 


KF372365 


KF372129 


KF372288 


KF372051 


DB8958 


13 


Guadarrama (Spain) 


40.71 N 


-4.14 E 


KF372207 


KF372366 


KF372130 


KF372289 


KF372052 


DB8959 


13 


Guadarrama (Spain) 


40.71 N 


-4.14 E 


KF372208 


KF372367 


- 


- 


- 


DB13472 


14 


Gudar (Spain) 


40.37 N 


-0.67 E 


- 


- 


KF372131 


- 


KF372053 


DB13473 


14 


Gudar (Spain) 


40.37 N 


-0.67 E 


KF372209 


KF372368 


KF372132 


KF372290 


KF372054 


DB13477 


15 


Penyagolosa (Spain) 


40.23 N 


-0.35 E 


KF372210 


KF372369 


KF372133 


KF372291 


KF372055 


DB13478 


15 


Penyagolosa (Spain) 


40.23 N 


-0.35 E 


KF372211 


KF372370 


KF372134 


KF372292 


KF372056 


DB 14060 


16 


Montseny (Spain) 


41.77 N 


2.44 E 


KF372212 


KF372371 


KF372135 


KF372293 


KF372057 


DB 14063 


16 


Montseny (Spain) 


41.77 N 


2.44 E 


KF372213 


KF372372 


KF372136 


KF372294 


KF372058 


DB 14064 


16 


Montseny (Spain) 


41.77 N 


2.44 E 


KF372214 


KF372373 


KF372137 


KF372295 


KF372059 


DB1759 


17 


Les Salines (Spain) 


42.42 N 


2.75 E 


KF372215 


- 


- 


- 


- 


DB5030 


18 


Meranges (Spain) 


42.43 N 


1.78 E 


KF372216 


KF372374 


KF372138 


KF372296 


KF372060 


DB11195 


18 


Meranges (Spain) 


42.43 N 


1.78 E 


KF372217 


KF372375 


KF372139 


KF372297 


- 


DB13718* 


19 


Luzenac (France) 


42.75 N 


1.75 E 


KF372218 


KF372376 


- 


- 


- 


DB 13461 


20 


Massif des Maures (France) 


43.24 N 


6.38 E 


KF372219 


KF372377 


KF372140 


KF372298 


KF372061 


DB 13460 


21 


Valle de Gilly (France) 


43.28 N 


6.46 E 


KF372220 


KF372378 


KF372141 


KF372299 


KF372062 


DB13430 


22 


Massif des Maures (France) 


43.38 N 


6.62 E 


KF372221 


KF372379 


KF372142 


KF372300 


KF372063 


DPMI 


23 


Viozene (Italy) 


44.14 N 


7.78 E 


KF372222 




KF372143 


KF372301 


KF372064 


DPM3 


23 


Viozene (Italy) 


44.14 N 


7.78 E 


KF372223 


KF372380 


KF372144 


KF372302 


KF372065 


DB15936 


24 


Monte Verita (Switzerland) 


46.16 N 


8.72 E 


KF372224 










DB16837 


25 


Bianzano (Italy) 


45.77 N 


9.92 E 


KF372225 


KF372381 


KF372145 


KF372303 


KF372066 


DB16838 


25 


Bianzano (Italy) 


45.77 N 


9.92 E 


KF372226 


KF372382 


KF372146 


KF372304 




DPM36 


26 


Peschiera del Garda (Italy) 


45.44 N 


10.67 E 


KF372227 


KF372383 


KF372147 


KF372305 


KF372067 


DPM37 


26 


Peschiera del Garda (Italy) 


45.44 N 


10.67 E 


KF372228 


KF372384 




KF372306 


KF372068 


DPM38 


26 


Peschiera del Garda (Italy) 


45.44 N 


10.67 E 


KF372229 


KF372385 


KF372148 


KF372307 


KF372069 
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Table 1 Geographical location and codes for the studied specimens of Podarcis muralis and the outgroup species 
P. liolepis (Continued) 



DB1399 


27 


Ostia Antica (Italy) 


41.75 N 


12.31 E 


KF372230 


KF372386 


KF372149 


KF372308 


KF372070 


DB5938 


28 


Paganico (Italy) 


42.19 N 


13.00 E 


KF372231 


KF372387 


KF372150 


KF372309 


KF372071 


DPM39 


29 


Majelletta (Italy) 


42.16 N 


14.13 E 


KF372232 


KF372388 


KF372151 


KF372310 


KF372072 


DPM40 


29 


Majelletta (Italy) 


42.16 N 


14.13 E 


KF372233 


KF372389 


KF372152 


KF37231 1 


KF372073 


DPM41 


29 


Majelletta (Italy) 


42.16 N 


14.13 E 


KF372234 








KF372074 


DPM15 


30 


Bassiano (Italy) 


41.55 N 


13.05 E 


KF372235 


KF372390 


KF372153 


KF372312 


KF372075 


DPM16 


30 


Bassiano (Italy) 


41.55 N 


13.05 E 


KF372236 


KF372391 


KF372154 


KF372313 


KF372076 


DPM18 


31 


Fondi (Italy) 


41.37 N 


13.33 E 


KF372237 


KF372392 


KF372155 


KF372314 


KF372077 


DPM19 


31 


Fondi (Italy) 


41.37 N 


13.33 E 


KF372238 


KF372393 


KF372156 


KF372315 


KF372078 


DPM25 


32 


Pollino National Park (Italy) 


39.93 N 


16.17 E 


KF372239 


KF372394 


KF372157 


KF372316 


KF372079 


DB5937 


33 


Pollino National Park (Italy) 


39.90 N 


16.19 E 


KF372240 


KF372395 


KF372158 


KF372317 


KF372080 


DPM21 


33 


Pollino National Park (Italy) 


39.90 N 


16.19 E 


KF372241 


KF372396 


KF372159 


KF372318 


KF372081 


DPM7 


34 


Fago del Soldato (Italy) 


39.35 N 


16.55 E 


KF372242 


KF372397 


KF372160 


KF372319 


KF372082 


DPM8 


34 


Fago del Soldato (Italy) 


39.35 N 


16.55 E 


KF372243 


KF372398 


KF372161 


KF372320 


KF372083 


DPM6 


34 


Fago del Soldato (Italy) 


39.35 N 


16.55 E 


KF372244 


KF372399 


KF372162 


KF372321 


KF372084 


DB13929 


35 


Ribnica (Slovenia) 


45.75 N 


14.77 E 


KF372245 


KF372400 


KF372163 


KF372322 


KF372085 


DB13930 


35 


Ribnica (Slovenia) 


45.75 N 


14.77 E 


KF372246 


KF372401 


KF372164 


KF372323 


KF372086 


DB13943 


36 


Donja Lamana Draga (Slovenia) 


45.51 N 


14.96 E 


KF372247 


KF372402 


KF372165 


KF372324 


KF372087 


DB13945 


36 


Donja Lamana Draga (Slovenia) 


45.51 N 


14.96 E 


KF372248 


KF372403 


KF372166 


KF372325 


KF372088 


DB15991 


37 


Zuce (Serbia) 


44.68 N 


20.55 E 


KF372249 


KF372404 


KF372167 


KF372326 


KF372089 


DB15992 


37 


Zuce (Serbia) 


44.68 N 


20.55 E 


KF372250 


KF372405 
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DB15993 
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44.68 N 


20.55 E 


KF372251 
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KF372169 


KF372328 


KF372091 


DB13309 
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Ci/iheBO (Serbia) 


43.34 N 


22.08 E 


KF372252 


KF372407 


KF372170 


KF372329 


KF372092 


DB16121 
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Road to Crnovska River (Serbia) 


42.38 N 


22.05 E 


KF372253 


KF372408 


KF372171 


KF372330 


KF372093 


DB16122 
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Road to Crnovska River (Serbia) 


42.38 N 


22.05 E 


KF372254 
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Velouxi (Greece) 


38.85 N 


21.38 E 


KF372258 


KF372413 
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Platanitsa, Achaia (Greece) 


37.97 N 


21.88 E 


KF372259 


KF372414 


KF372176 


KF372336 


KF372099 


DPM35 


44 


Platanitsa, Achaia (Greece) 


37.96 N 


21.91 E 


KF372260 




KF372177 


KF372337 


KF372100 


DB16936 


45 


Mainalo, Pelloponnisos (Greece) 


37.39 N 


22.46 E 


KF372261 


KF372415 








DB16933 


46 


Kisavos mt (Greece) 


39.63 N 


22.64 E 


KF372262 


KF372416 








DB16938 


47 


Samothraki isl. (Greece) 


40.45 N 


25.59 E 


KF372263 


KF372417 


KF372178 


KF372338 




DB16310 


48 


Kapakli (Turkey) 


41.89 N 


27.35 E 


KF372264 


KF372418 


KF372179 


KF372339 


KF372101 


DB16307 


49 


Kapakli (Turkey) 


41.91 N 


27.36 E 
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50 


Derekoy (Turkey) 


41.96 N 


27.39 E 


KF372267 


KF372421 


KF372182 


KF372342 


KF372104 


DB 16482 


50 


Derekoy (Turkey) 


41.96 N 


27.39 E 




KF372422 


KF372183 


KF372343 


KF372105 


DB16322 


51 


Derekoy (Turkey) 


41.97 N 


27.42 E 


KF372268 


KF372423 


KF372184 


KF372344 


KF372106 


DB16332 


51 


Derekoy (Turkey) 


41.97 N 


27.42 E 


KF372269 


KF372424 


KF372185 


KF372345 


KF372107 


DB 16485 


51 


Derekoy (Turkey) 


41.97 N 


27.42 E 


KF372270 


KF372425 


KF372186 


KF372346 


KF372108 


DB16318 


52 


Derekoy (Turkey) 


41.96 N 


27.45 E 




KF372426 


KF372187 


KF372347 


KF372109 
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Table 1 Geographical location and codes for the studied specimens of Podarcis muralis and the outgroup species 
P. liolepis (Continued) 

DB16341 53 Pi narozii Turkey) 41.77 N 34.04 E KF372271 KF372427 KF372188 KF372348 KF372110 

DB16340 53 Pinarozu (Turkey) 41.77 N 34.04 E - KF372428 KF372189 KF372349 KF3721 1 1 

DB16451 53 Pinarozu Turkey) 41.77 N 34.04 E KF372272 KF372429 KF372190 - KF372112 

Specimen and code locality along with geographic coordinates are reported. Genbank Accession Numbers for each gene are provided. *outgroup species: 
P. liolepis. 



FJ460597). Haplotype phase of nuclear genes was deter- 
mined using PHASE version 2.1.1. [35]. Estimations were 
carried out under a model with recombination (-MR0 
option), with the initial 1000 iterations discarded as 
burn-in, 1 as thinning interval and 1000 post-burnin it- 
erations. Three independent runs were conducted for 
each gene. After monitoring the goodness of fit for each 
run according to the programs manual, we accepted 
haplotype reconstructions which both i) yielded the 
same result in each of the three runs, and ii) had associ- 
ated phase posterior probabilities of at least 0.75 in the 
average of the three runs. For each nuclear gene, the 
possible occurrence of recombination events was assessed 
using the Pairwise Homoplasy Index (PHI) test [36] 
implemented in SPLITSTREE 4.11 [37]. 

In order to assess patterns of genetic diversity and 
evaluate whether southern peninsulas harboured most of 
the diversity of P. muralis, we estimated the overall gen- 
etic diversity of the species and its distribution across 
populations and regions. For each gene we computed the 
following summary statistics of genetic diversity: number 
of segregating sites S, nucleotide diversity jt, number of 
haplotypes H, and haplotype diversity Hd, both overall 
and for specific groups or clades defined within the spe- 
cies (see Results and Discussion) with DNASP 5 [38]. 
Comparisons of diversity measures can be greatly affected 
by different sample sizes among groups. In order to ac- 
count for this, we followed a resampling approach: given 
the lowest sample size (except 0 or 1) amongst the various 
groups being compared (s), for each group with sample 
size > s we resampled 100 sets of s sequences, calculated 
diversity measures for each of these 100 sets and took the 
average for each diversity measure. This procedure allows 
a straightforward comparison of diversity measures de- 
rived from groups or clades with distinct sample sizes. 
Resampling was performed with the aid of a series of 
scripts written in Python 2.7.1 (available from the authors 
upon request) and taking advantage of DNASP "batch 
mode" calculations. Average uncorrected genetic distances 
between mitochondrial clades were assessed using MEGA 
5 [39]. Furtehrmore, to assess genetic differentiation 
between clades we calculated F ST (based on ^-distance) 
and its significance by performing 1000 permutations in 
ARLEQUIN 3.5.1.2 [40]. 

Phylogenetic relationships among mitochondrial hap- 
lotypes were inferred using Maximum Likelihood (ML) 



and Bayesian (BA) methods. We tested three outgroups, 
P. liolepis, P. bocagei, and P. sicula (both combined and 
separately) through preliminary ML and BA tree searches. 
The outgroup effect on the ML and BA topologies was 
limited and differences between trees calculated using dif- 
ferent outgroups were limited to unsupported nodes (not 
shown). Thus, we present and discuss phylogenetic ana- 
lyses using P. liolepis as an outgroup (Genbank Accession 
numbers for P. bocagei cytblnd4: DQ081139/DQ081153; 
for P. sicula cytb/nd4\ KF372034/KF372035; for P. liolepis 
see Table 1). The ML search and the model selection were 
carried out in TREEFINDER version October 2008 [41]. 
We selected the best partitioned model of nucleotide 
substitution under the corrected Akaike s Information Cri- 
terion (AICc) by comparing the AICc scores of seven dif- 
ferent partition schemes (under a fixed ML topology): (i) 
by gene fragments (cytb /nd4+tRN As), (ii) by genes {cytbl 
nd4/tRNAs); (iii) by coding regions (cytb+nd4/tRN As), 
(iv-v) by codon position alone (two schemes: lst/2nd/3rd 
and lst+2nd/3rd), and (vi-vii) by codon position, in com- 
bination with gene partitions (for a total of seven or five 
partitions, depending on the codon partitioning scheme). 
The optimal model had seven partitions, three codon par- 
titions for each coding gene and one partition for tRNAs 
{cytb: HKY for the 1st, the 2nd, and the 3rd codon posi- 
tions; nd4: TN+ G for the the 1st, and HKY+G for the 
2nd and the 3rd codon positions; tRNAs: HKY+G, each 
model with frequency and rate heterogeneity parameters 
optimisation). We performed a global tree search using 
100 random start trees generated through equidistant 
random walks of random nearest-neighbour-interchanges 
(NNI) starting from the centre tree obtained by a simple 
ML search. Support for the nodes (BP) was evaluated with 
1000 bootstrap replicates. Additionally, in order to evalu- 
ate the total mitochondrial variation known for this spe- 
cies, we carried out a ML analysis following the same 
procedure as above (the best model in this case was the 
HKY with a single partition), using a dataset including 
both the cytb sequences generated in this study (N=82) 
and those available in GenBank for native populations of 
P. muralis (N=103). This dataset provides even greater 
coverage of the whole species range (see Figure 1). The 
obtained phylogenetic trees were visualised and edited 
using FigTree v. 1.3.1 (available at http://tree.bio.ed.ac. 
uk/software/figtree/). The geographic pattern of haplo- 
group distribution and the relationships between 



Salvi et al. BMC Evolutionary Biology 2013, 13:147 
http://www.biomedcentral.com/1471 -21 48/1 3/1 47 



Page 8 of 18 



haplotypes help to disentangle the two alternative sce- 
narios predicted by the southern refugia model or by a 
model of Mediterranean and extra-Mediterranean refugia. 
For this purpose, we mapped the distribution of 
haplogroups identified by the phylogenetic trees across 
our range-wide sampling. 

Bayesian analyses were carried out using the BEAST 
v. 1.7.4 package [42] for two purposes: i) to estimate a phyl- 
ogeny of the combined mitochondrial DNA dataset based 
on Bayesian inference; ii) to obtain estimates of the 
mtDNA time to the most recent common ancestor 
(TMRCA) of P. muralis and groups detected within the 
species. The choice of an appropriate tree prior for this 
dataset is somewhat problematic: on one hand, the bulk of 
our data is intraspecific, hence a "speciation" prior such as 
the Yule process may not be adequate. On the other hand, 
because the presence of an outgroup is mandatory and due 
to the high levels of geographic substructure detected in 
our focal species (see Results), our data are not totally 
amenable to analyses under a basic coalescent process. 
Therefore, we performed these analyses under both a Yule 
process and a coalescent process with constant size, to 
evaluate the effect of prior selection on our estimates of 
phylogeny and TMRCAs. To avoid overparametrization, 
we implemented a model with a simpler partition scheme, 
defined by gene fragments (two-partitions), and selected 
under AICc in TREEFINDER. The best evolutionary model 
for our dataset assumed the HKY+G and the TIM2+G+I 
substitution models for the cytb and nd4 fragments, re- 
spectively. The TIM2 model was manually specified within 
the BEAUti input file [42]. To infer mtDNA TMRCAs, we 
assumed a relaxed molecular clock for both cytb and nd4, 
applying an uncorrected lognormal model [43]. We used 
the available rates proposed for Podarcis for the longer 
fragment, nd4, while providing a large uninformative prior 
for the cytb substitution rate (parameter ucld.mean). In the 
case of nd4 we assumed that differences accumulate at a 
mean rate of 2.26% per million years, obtained as the aver- 
age between divergence rates proposed for Podarcis for this 
fragment (ranging from 1.74% to 2.78% per million years; 
[44]). In order to incorporate this uncertainty into TMRCA 
estimates, we defined a normal prior distribution on the 
mutation rate with mean 1.13%/MY and with 95% of the 
probability distribution encompassed between the two 
aforementioned limits (0.87 and 1.39%, respectively). For 
each tree prior, BEAST was run twice, with 100 million iter- 
ations per run, sampling every 10000 steps. TRACER v. 1.5 
(available at http://beast.bio.ed.ac.uk/Tracer) was used to 
visualize the results and assess if the effective sample size 
of estimated parameters was satisfactory. We discarded 
the initial 10% of sampled trees as burn-in. Upon confirm- 
ation of convergence, the two outputs for each tree prior 
were combined and subsequent data were retrieved from 
the combined data set. 



For the following analyses, we considered as clades 
groups of haplotypes which are reciprocally monophyletic 
at the mitochondrial loci with the ML and BA tree recon- 
structions; displaying geographical contiguity; and diver- 
gent by more than 1% with its closer relative in line with 
previous studies ([24,27,28]; see also [45]). While we ac- 
knowledge that such criteria could be considered arbitrary, 
to further test that these groups correspond to evolution- 
ary significant lineages, ESU sensu [46], we confirmed they 
show significant divergence of allele frequencies at nuclear 
loci as assessed by Fst statistics (see above). 

To assess how genetic variance at the combined cytb I 
nd4 sequences was hierarchically distributed among and 
within mitochondrial clades, we performed a spatial ana- 
lysis of molecular variance (SAMOVA) with SAMOVA 
1.0 [47]. We pooled samples from neighbouring local- 
ities in accordance with the mitochondrial clades previ- 
ously identified. Geographic coordinates for each clade 
were calculated as the geographic centroid of member 
localities. Clades represented by less than two individ- 
uals were not included in the analysis. SAMOVA was 
run for 10,000 iterations from each of 100 random 
initial conditions, and testing the predefined number of 
groups (K) from 2 to 14. 

Phylogenetic relationships among nuclear haplotypes 
were inferred by median joining networks [48] using the 
software NETWORK 4.6.1.0 (available at http://www. 
fluxus-engineering.com/sharenet.htm). 

Finally, we specifically tested the hypothesis of a recent 
expansion of P. muralis from southern glacial refugia 
against the alternative scenario of Pleistocene survival of 
the species in both southern and northern refugia. The 
first hypothesis implies that northern localities belong to 
the same genetic pool as the southern region from which 
they originated. As such, sequence divergence between 
northern and southern localities will be, in most cases, 
significantly lower than the divergence observed between 
clades and possibly similar to intraclade divergence. By 
contrast, the latter hypothesis implies that northern re- 
gions have hosted surviving populations which have 
evolved independently from southern ones. Thus, we ex- 
pect the mean sequence divergence between northern 
and southern localities to be similar to mean divergence 
among clades. This rationale provides a statistical frame- 
work to explicitly compare the two competing hypoth- 
eses. A requisite of this test is that mean between-clade 
sequence divergence is significantly higher than within- 
clade divergence. Therefore, as a first step, we performed 
a permutation test to evaluate if between-clade ^-dis- 
tances (D bc ) were significantly higher than within-clade 
distances (D wc ). For this purpose, we randomly shuffled 
individuals among clades 1000 times in order to obtain 
the empirical distribution of the difference D bc - D wc . 
The corresponding ^-value was then obtained as the 
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percentage of times for which the observed D bc - D wc 
was smaller than random. This test was performed for 
the combined mitochondrial DNA data set, as well as 
for the three nuclear loci in separate, always using the 
units defined on the basis of the mtDNA phylogeny. In 
continuation, for each northern locality and for each 
locus, we evaluated if the sequence divergence from 
southern localities fell within the "between-clade" or the 
"within-clade" level of differentiation. We classified lo- 
calities as "northern" or "southern" according to biogeo- 
graphic barriers separating central Europe from major 
glacial refugia in southern peninsulas as defined by pre- 
vious literature [1,4,5,7,15]: localities north or east to the 
Pyrenees (locality 1), north to the Appennine chain (lo- 
calities 25-26), west or east to the Alps (localities 20-23 
and 35-36, respectively), and in the central Balkans (lo- 
calities 37-39) were therefore considered as "northern". 
For each northern locality we performed two different 
tests: first, we examined its mean ^-distance from all 
southern localities (localities 2-16, 18, 27-34, and 40- 
46). Second, to account for the possibility that northern 
localities may have been preferably colonized from adja- 
cent refugia, we also examined the mean ^-distance of 
each northern locality from contiguous southern regions 
only (locality 1 vs localities 2-16, and 18; localities 20- 
26 vs localities 27-34; localities 37-39 vs. localities 40- 
46; localities 35-36 were compared either with localities 
27-34 or with localities 40-46, given their intermediate 
geographic location). These comparisons were performed 
according to the following permutation procedure: indi- 
viduals were shuffled among localities at random across 
1000 permutation cycles. In each permutation, we calcu- 
lated the mean ^-distance between each northern locality 
and the set of southern localities being compared (D b 
where i is the code of the northern locality), the mean 
between-clade distance (D bc ) and the mean within-clade 
distance (D wc ). For each northern locality z, we then 
assessed whether the difference between D bc and D t (D bc - 
D t ) was equal or higher than the distance we observed in 
the real data. The same was done for the difference be- 
tween D t and D wc {D t - D wc ). The proportion of permuta- 
tions verifying each of these conditions was used as a 
j?-value to assess if the observed differences were signifi- 
cant or could be due to chance alone. 

Results 

We obtained 82 cytb mitochondrial sequences of 411 
base pairs (bp), 79 nd4 mitochondrial sequences of 873 
bp (75 individuals sequenced for both fragments), 156 
sequences for the nuclear gene fragment mclr, and 154 
for acm4 and pdc with 688bp, 423 bp, and 342 bp, re- 
spectively. GenBank accession numbers for the sequences 
generated in this study are reported in Table 1. Multiple 
sequence alignment did not require gap positions in any 



of the studied genes. The phi test did not find statistically 
significant evidence for recombination within nuclear 
gene fragments (p > 0.05). Sample size and localities, and 
clade assignment are reported in Table 2. 

The phylogenetic relationships inferred from the ML 
tree based on the combined mitochondrial sequences 
showed 17 well supported clades (BP>90) with a geo- 
graphic coherence (Figure 1A). The Bayesian tree showed 
an identical topology at supported nodes (not shown) with 
the same 17 clades supported with a posterior probability 
of 1.00 (Figure 1A), independently of the choice of tree 
prior. Clade 1 grouped haplotypes from Northwest France 
and Southeast Pyrenees (Spain). Haplotypes included in 
clades 2 and 3 occurred in populations from Northern 
and Central Spain, respectively. Clade 6 is restricted to 
Provence, France. Five different clades occurred in Italy: 
clade 5 in Northern Italy beneath the Alpine arch; clades 9 
and 10 in central Italy; and clades 7 and 8 in the Calabrian 
Peninsula. Slovenian populations clustered in clade 11 and 
samples from Serbia in clade 4. In the Balkan Peninsula, 
clades 12, 13, and 14 were distributed in Thessaly, Central 
Greece, and Peloponnese regions, respectively. Clade 16 is 
restricted to Samothraki Island and clades 15 and 17 are 
distributed in Turkish Thrace and North- Western Anato- 
lia, respectively. Clades 12 and 16 were formed by single 
haplotypes, and as such they cannot be strictly considered 
as clades. However, we use the term clade' for simplicity, 
given their divergence from other clades, and, in the case 
of clade 12, also because this is indeed formed by multiple 
haplotypes in the analyses including Genbank sequences 
(see below and Figure IB). The average uncorrected gen- 
etic distance between mitochondrial clades was 4.5 and 
3.1% for cytb and nd4 respectively (ranging from 1 to 7.4% 
for cytb and 0.9 to 4.7% for nd4) and both loci show sig- 
nificant F ST values (0.86 and 0.83 for cytb and nd4 respect- 
ively; P=0.001). Clades 15, 16, and 17 showed the lowest 
pairwise distance values at both cytb (1%) and nd4 (0.9%), 
while the minimum value of pairwise genetic distances 
across other clades was 1.8% for both genes. Phylogenetic 
relationships between Iberian clades West of the Pyrenees 
(clades 2 and 3), between southern Italian clades (clades 7 
and 8), between Slovenian and Balkan clades (clades 11 
and 12), and among clades from Turkey and Samothraki 
(clades 15-17) were well supported in both ML and BA 
trees, while the relationships between clades 1-6 was sup- 
ported in BA but not in ML analyses (BP<60). The rela- 
tionships between the remaining clades were not resolved. 

The ML analysis combining the 84 cytb sequences 
generated in this study and 103 sequences published 
from previous studies showed 23 well supported clades 
(Figure IB). These included the 17 clades identified 
based on the combined cytblnd4 ML tree and are num- 
bered according to the combined ML tree (Figure 1A) 
and named according to previous studies [22,26,27]. 



Table 2 Summary statistics of genetic diversity for each gene, clade, and region 
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3 


5 


0.93 


4.33 


6 


342 


3 


3 


0.60 


4.05 








754 


2.29 


2 


1.00 


3.04 




391 


0.62 


1.62 


0.62 


1.64 




576 


0.29 


1.29 


0.29 


0.50 




323 


1.44 


1.97 


0.97 


4.46 




248 


1.92 


2.29 


0.58 


3.89 


8 


34 


3 


873 


11 


3 


1.00 


9.73 


3 


411 


10 


2 


0.67 


17.64 


6 


683 


1 


2 


0.33 


0.58 


6 


423 


2 


3 


0.80 


3.30 


6 


342 


1 


2 


0.33 


1.35 








754 


7.54 


2 


1.00 


10.00 




391 


7.5 


1.75 


0.75 


19.84 




576 


0.29 


1.29 


0.29 


0.50 




323 


1.05 


1.78 


0.78 


3.25 




248 


0.62 


1.62 


0.31 


1.25 


9 


27, 28, 30, 31 


6 


873 


18 


6 


1.00 


10.52 


6 


411 


9 


3 


0.73 


12.17 


12 


683 


4 


5 


0.58 


1.16 


12 


423 


1 


2 


0.17 


0.52 


12 


342 


2 


3 


0.67 


3.44 








754 


8.3 


2 


1.00 


11.01 




391 


4.94 


1.78 


0.78 


13.07 




576 


0.63 


1.54 


0.54 


1.09 




323 


0.21 


1.21 


0.21 


0.65 




248 


1.45 


2.33 


0.62 


3.25 


10 


29 


2 


873 


2 


2 


1.00 


2.65 


3 


411 


1 


2 


0.67 


1.76 


4 


683 


2 


2 


0.50 


1.74 


4 


423 


1 


2 


0.50 


1.55 


6 


342 


1 


2 


0.53 


2.16 








754 












391 


0.67 


1.67 


0.67 


1.77 




576 


1 


1.5 


0.50 


1.74 




323 


0.49 


1.49 


0.49 


1.52 




248 


0.91 


1.91 


0.52 


2.11 


11 


35, 36 


4 


873 


6 


4 


1.00 


4.42 


4 


411 


1 


2 


0.50 


1.32 


8 


683 


2 


3 


0.71 


1.49 


8 


423 


0 


1 


0.00 


0.00 


8 


342 


1 


2 


0.43 


1.74 








754 


3.28 


2 


1.00 


4.35 




391 


0.52 


1.52 


0.52 


1.38 




576 


0.85 


1.71 


0.71 


1.48 




323 


0 


1 


0.00 


0.00 




248 


0.76 


1.76 


0.41 


1.67 


12 


46 


1 


873 
754 


- 


- 


- 


- 


1 


411 
391 


- 
- 


- 
- 


- 
- 


- 
- 


0 


- 
- 


- 
- 


- 
- 


- 
- 


- 
- 


0 


- 
- 


- 
- 


- 
- 


- 
- 


- 
- 


0 


- 
- 


- 
- 


- 
- 


- 
- 


- 
- 


13 


40,41 


2 


873 
754 


3 


2 


1.00 


3.98 


2 


411 
391 


2 


2 


1.00 


5.29 


4 


683 
576 


1 

0.52 


2 
1.52 


0.50 
0.52 


0.87 
0.90 


4 


423 
323 


2 
1.05 


3 
1.88 


0.83 
0.88 


3.10 
3.25 


4 


342 
248 


0 


1 


0.00 


0.00 


14 


42-45 


3 


873 


16 


3 


1.00 


14.15 


4 


411 


10 


4 


1.00 


14.11 


4 


683 


1 


2 


0.50 


0.87 


6 


423 


0 


1 


0.00 


0.00 


6 


342 


0 


1 


0.00 


0.00 








754 


10.17 


2 


1.00 


13.49 




391 


5.26 


2 


1.00 


13.92 




576 


0.51 


1.51 


0.51 


0.89 




323 


0 


1 


0.00 


0.00 




248 


0 


1 


0.00 


0.00 


15 


48-52 


9 


873 


11 


9 


1.00 


4.86 


7 


411 


0 


1 


0.00 


0.00 


18 


683 


0 


1 


0.00 


0.00 


18 


423 


2 


3 


0.45 


2.45 


18 


342 


0 


1 


0.00 


0.00 








754 


3.85 


2 


1.00 


5.11 




391 


0 


1 


0.00 


0.00 




576 


0 


1 


0.00 


0.00 




323 


0.8 


1.45 


0.45 


2.48 




248 


0 


1 


0.00 


0.00 



Table 2 Summary statistics of genetic diversity for each gene, clade, and region (Continued) 



16 


47 


1 


873 










1 


411 










2 


683 


0 


1 


0.00 


0.00 


2 


423 


0 


1 


0.00 


0.00 


0 


















754 












391 


_ 


_ 


_ 


_ 




576 


_ 


_ 


_ 


_ 




323 


_ 


_ 


_ 


_ 




_ 


_ 


_ 


_ 


_ 


17 


53 


3 


873 


2 


3 


1.00 


1.77 


2 


411 


0 


1 


0.00 


0.00 


6 


683 


0 


1 


0.00 


0.00 


4 


423 


2 


2 


0.67 


4.13 


6 


342 


0 


1 


0.00 


0.00 








754 


1.26 


2 


1.00 


1.67 




391 












576 


0 


1 


0.00 


0.00 




323 


1.34 


1.67 


0.67 


4.15 




248 


0 


1 


0.00 


0.00 






79 


873 


142 


74 


1.00 


28.21 


81 


411 


83 


42 


0.97 


40.05 


156 


688 


24 


21 


0.37 


0.92 


154 


423 


11 


13 


0.57 


2.40 


154 


342 


13 


12 


0.73 


6.67 


Region 


Clades 






























































Balkans 


12, 13, 14 


7 


873 


53 


7 


1 


27.92 


8 


411 


38 


8 


1 


36.66 


10 


683 


2 


3 


0.38 


0.69 


12 


423 


3 


4 


0.56 


1.97 


10 


342 


0 


1 


0 


0.00 








754 












391 












576 












323 












248 










Italy 


7, 8, 9, 10 


14 


873 


60 


14 


1 


27.49 


16 


411 


43 


10 


0.94 


35.45 


28 


683 


8 


8 


0.44 


0.99 


28 


423 


3 


5 


0.57 


2.19 


30 


342 


6 


6 


0.62 


3.24 








754 


47.88 


7 


1 


27.25 




391 


34.87 


6.48 


9.42 


35.35 




576 


2.6 


3.31 


0.41 


0.90 




323 


2.58 


3.85 


0.56 


2.17 




248 


2.76 


3.47 


0.62 


3.20 


Iberia 


1*, 2, 3 


24 


873 


32 


22 


0.99 


15.29 


25 


411 


21 


9 


9.78 


18.94 


50 


683 


8 


6 


0.19 


0.56 


46 


423 


2 


3 


0.13 


0.40 


46 


342 


3 


3 


0.58 


4.56 








754 


24.84 


6.89 


0.99 


14.75 




391 


17.66 


4.7 


0.79 


19.25 




576 


1.61 


2.06 


0.20 


0.56 




323 


0.77 


1.77 


0.14 


0.45 




248 


2.68 


2.69 


0.58 


4.62 



N, sample size (in number of gene copies), ns number of sites, S, number of segregating sites, h number of haplotypes, Hd haplotype diversity, n nucleotide diversity. For the diversity measures S, h, Hd, and n values 
calculated under a resampling approach are also reported below values calculated on the original datasets. For the number of sites, for each case the first row represents the size of the alignment while the second 
row represents the number of sites after discarding Ns or unphased positions; "-" means that no resampling was conducted for that set. ^excluding locality 1. 
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Previous studies had identified 12 of these mitochondrial 
mtDNA lineages - eight across Italy ("lAa", "lAc", "lAd", 
"lAe", "IB", "2a", "2b", and "2c"; [26]), one from the island 
of Elba [27], and one from each of Western France, 
Eastern France, and the Balkans [23]. We separated the 
clades "1A" and "IB" listed in [26] into their constituent 
sub-clades, since clade "1A" was not supported by phylo- 
genetic analyses and clade "IB" from southern Italy is 
composed of distinct units (from Calabria, Apulia and 
Campania), which are as differentiated as the other clades. 
In our assessment, ten previously undetected clades were 
identified (clades 2,3, 8, and 11-17 in Figure 1A). Add- 
itionally, a new lineage was identified in East Macedonia 
(clade 18) that had not been previously described as a 
group despite being comprised of sequences retrieved 
from GenBank. 

MtDNA TMRCA estimates within P. muralis were 
roughly similar independently of the choice of tree prior 
in BEAST analyses (Additional file 3). The TMRCA of P. 
muralis was estimated at 2.47-2.67 million years ago, 
thus coinciding with the Plio- Pleistocene boundary. The 
TMRCA estimates of the various clades bound the 
period from Early Pleistocene to Late Pleistocene. 

SAMOVA analyses showed a lack of genetic structure 
above the clade level (Additional file 4). No clade group- 
ing explained the genetic variation observed in P. 
muralis significantly better than other grouping options, 
as indicated by the uniform increase of F C t values to- 
ward higher values of K. 

The phylogenetic networks depicting the relationships 
between acm4, mclr, and pdc haplo types are shown in 
Figure 2. The acm4 and mclr networks showed a pattern 
of low divergence and a lack of apparent phylogeographic 
structure, although both show significant F ST values (0.37 
and 0.24 respectively; P=0.001) between mtDNA-defined 
clades. In both loci the most common haplotypes were 
found in most localities (with the exception of localities 47 
and 53 for acm4, while locality 46 is not represented in 
the dataset of these two genes, see Table 1 and 2) and de- 
rived haplotypes divergent by one or two mutational steps 
had a more restricted distribution. On the other hand, evi- 
dence of geographic association of nuclear haplotypes was 
more prominent for the pdc locus (Figure 2), leading to a 
higher and also significant Fst value (F S t = 0.56; P=0.001). 
The structure of the network showed two most frequent 
haplotypes, each exhibiting derived variants: the haplo- 
group placed in the left side of the network was restricted 
to the Balkans, Turkey, Slovenia, and Provence (France), 
while the haplogroup on the right side included all the 
haplotypes from Italy. Haplotypes sampled in the 
Iberian Peninsula and Eastern France belonged to both 
haplogroups. 

Estimates of genetic diversity «S, n, H, and Hd for each 
gene, clade, and geographic region are given in Table 2. 




Figure 2 Phylogenetic relationships among nuclear haplotypes. 

Haplotype median joining networks based on sequences of the 

nuclear genes ocm4, mclr, and pdc. Haplotypes are represented by 

circles with area proportional to their frequency and coloured 

according to mitochondrial clades as defined by the Maximum 

Likelihood tree (see Figure 1). 
k ) 

Nucleotide and haplotype diversity of mitochondrial 
genes was higher in Northern, Central, and Southern 
Italy (clades 5, 9, and 8, respectively) and in clade 14 
from Southern mainland Greece and the Peloponnese. 
At the nuclear gene level, an overall pattern of low gen- 
etic diversity was found. Higher diversity values were 
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observed in populations from Northern Italy (clade 5), 
whereas no clear trend was apparent across putative 
southern refugial areas. 

Permutation tests showed that between-clade ^-dis- 
tance is significantly higher than within-clade distances 
for all loci (D bc /D wc for mitochondrial DNA: 0.0347/ 
0.0052; for acm4: 0.0051/0.0022; for mclr: 0.0029/ 
0.0013; for pdc: 0.0097/0.0043; P=0.001; see Additional 
file 5). This allowed us to perform tests aiming at evalu- 
ating whether the genetic distance between northern 
and southern localities is consistent with the distance 
typically observed either between or within clades. Re- 
sults from these analyses showed that mean pairwise dis- 
tances between northern and southern localities are 
significantly different from within-clade differentiation, 
but not significantly different from between-clade dis- 
tance (except locality 1), independently of the locus used 
for these comparisons (for more details see Additional 
file 5). Similar results were obtained when comparing 
northern localities to the sub-set of contiguous southern 
locations, although in this case this pattern is strong at 
mitochondrial but less evident at nuclear loci. These 
results support the hypothesis that populations of P. 
muralis from northern regions correspond to distinct 
evolutionary units from those occurring in southern 
peninsulas, and therefore that they have likely originated 
in situ rather than as a result of postglacial expansions 
of southern populations. 

Discussion 

The common wall lizard Podarcis muralis exhibited a 
complex phylogeographic pattern with multiple divergent 
mtDNA clades across its range. In total 23 alio- parapatric 
clades were identified, most of them with a restricted 
distribution in the Iberian, Italian, and Balkan peninsulas, 
as well as in Western and Eastern France, Slovenia and 
Austria, Central Balkans, and North- Western Anatolia. 
The detection of previously unidentified lineages was evi- 
dently due to the extended sampling scheme of this study 
as compared to earlier assessments, which were either fo- 
cussed on a regional scale [27,28], or did not sample all 
clades endemic to the Iberian, Balkan or Anatolian penin- 
sulas [24]. Although the definition of clades is somewhat 
arbitrary, the lack of resolution of relationships between 
monophyletic haplogroups, their high divergence, their 
geographic coherence, and results from the SAMOVA 
analysis (Figure 1A,B; Additional file 4) indicate that a 
grouping into fewer lineages would not be practical and 
does not describe accurately the mitochondrial DNA di- 
versity exhibited by P. muralis. 

The observed phylogeographic structure of P. muralis 
does not match the current subspecific division of this 
species [22,23], as already reported by [27] and [28] for 
the 11 subspecies occurring in Italy. Regarding the 



subspecies found outside Italy, many divergent clades 
fall into P. m. brogniardii (clades 1-3, 5, and 6) and P. 
m. albanica (clades 12-14, and 18). Though a taxo- 
nomic revision is beyond the scope of this study, our re- 
sults indicate that the current intraspecific taxonomy of 
P. muralis reflects local phenotypic varieties rather than 
distinct evolutionary units. 

While mitochondrial data showed high variability and 
a clear geographic structure of genetic diversity, nuclear 
data showed both relatively low genetic variability and 
shallow divergence among populations. This finding is in 
agreement with allozyme data [49]; see also [26] from 
which low levels of polymorphism and genetic differenti- 
ation were identified (number of alleles per locus, A = 
1.07; proportion of polymorphic loci, P = 0.07; observed 
heterozygosity, H Q = 0.028; genetic distance D Nm = 0.036; 
average values from [49]) among 23 populations of P. 
muralis attributable to eight mitochondrial lineages de- 
scribed in our study (clades 1, 2, 3, 5, 9, 14, 22, and 23). 
Despite low differentiation in nuclear markers, F ST 
values between mitochondrial clades are significant for 
all three nuclear markers and the mean between-clade 
pairwise distance is significantly different from that ob- 
served within clades, indicating that they are significant 
evolutionary lineages sensu [46]. Furthermore, even though 
overall nuclear data showed low genetic variation, a weak 
geographic association of nuclear haplotypes is apparent in 
sequence data. At the acm4 and mclr loci, besides a com- 
mon widespread haplotype, some closely related derived 
haplotypes, relatively frequent, were found in the Balkans, 
Turkey, Slovenia, and Eastern France. This geographic 
association was more evident at the pdc locus, where we 
observed two main haplotype groups: one including all 
haplotypes from Italy, and the other comprising haplotypes 
from the Balkans, Turkey, Slovenia, and Eastern France. 
Both haplogroups also occurred in Spain and Western 
France. Such shallow phylogeographic patterns at nuclear 
loci suggest that the divergence between mitochondrial 
clades has not been long enough for them to reach recipro- 
cal monophyly at the nuclear level [50-52]. According to 
mitochondrial DNA TMRCA estimates, divergence among 
lineages is estimated to have occurred since the Early 
Pleistocene. Taking into account the slower evolutionary 
rates of nuclear genes, the genetic pattern observed at the 
mitochondrial and nuclear level is compatible with a sce- 
nario of Pleistocene divergence among lineages. Yet, the 
use of more variable markers such as microsatellites is 
needed for a full assessment of the nuclear pattern of vari- 
ation within this species. 

The finding of conspicuous genetic divergence among 
geographic clusters of populations has provided evidence 
for the existence of cryptic species within other Podarcis 
taxa including the P. erhardii [53], P. tiliguerta [31], and 
P. hispanica [54] species complexes. Given the extensive 
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genetic structure observed in P. muralis and the lack of 
evidence for syntopy between mitochondrial lineages 
(except in four populations in Italy, see Figure 1), it might 
be suggested that some form of barrier to gene exchange 
has existed between distinct lineages. Nevertheless, we 
have no evidence that P. muralis constitutes a species 
complex. Pairwise uncorrected ^-distances between clades 
ranges from 1.0 to 7.4% in cytb and from 0.9% to 4.7% in 
nd4, which is a level of divergence typically found within 
(and not between) recognised species of the genus (e.g. P. 
vaucheri, [44]). Further, there is no evidence for high dif- 
ferentiation in nuclear markers which would be indicative 
of long-term independent evolutionary trajectories. While 
a low divergence level among genetic lineages is not 
necessarily a synonym of absence of reproductive isola- 
tion, the available evidence suggest that, contrary to other 
Podarcis forms, P. muralis is a single species and not a 
species complex. 

This does not mean, however, that the genetic land- 
scape of the species is uniform. On the contrary, genetic 
variation in P. muralis is structured in many distinct 
lineages across its range. The phylogeographic pattern 
observed, with many reciprocally monophyletic allo- 
parapatric lineages having a Pleistocene divergence, sug- 
gests a scenario of long-term isolation in multiple ice-age 
refugia [1,55]. Therefore, the geographic distribution of 
these genetic lineages is likely to reflect the refugial struc- 
ture of this species and indicate putative areas which 
allowed the long-term maintenance of its genetic diversity. 
Assessments of many widespread European species have 
recovered a classic pattern of higher genetic diversity and 
multiple genetic lineages within southern European penin- 
sulas [2]. On one level the pattern of P. muralis is not an 
exception: more than three quarters of all lineages are en- 
demic to the Iberian (2 and 3), the Italian (7-10, 19-23) 
and the Balkan peninsulas (12-14, 18). The simplest ex- 
planation for this phylogeographic pattern is one of a large 
number of refugia occurring not only across southern 
peninsulas, as pioneer works in phylogeography have 
demonstrated at a continental scale (e.g. [1]), but also 
within each peninsula, conforming to the "refugia within 
refugia" paradigm [56]. Such compartmentalization is par- 
ticularly evident in Italy. Multiple, related lineages occur 
in South Italy, whereas other, unrelated lineages occur in 
central and northern regions. These latter lineages from 
central and northern regions were considered by [27] as 
belonging to a main lineage (clade "1A") leading to the in- 
ference of a main refugial area north of Campania and 
Apulia. However, though this clade was not supported in 
any of their phylogenetic analyses (BP < 56, BPP = 0.54), 
in our analyses its constituent lineages are statistically well 
supported, unrelated, well differentiated one from the 
other, and parapatrically distributed. Thus, a scenario of 
multiple independent refugia in Central and Northern 



Italy appears to fit the data better. A similar reasoning 
suggests the existence of multiple independent refugia in 
the southern Balkan Peninsula. On the other hand, the 
southern Italian Peninsula would have acted as an inde- 
pendent biogeographic compartment with four related 
lineages arising from distinct refugia within this main 
refugial area. The same situation is observed in the west- 
ernmost and easternmost areas of the species range - in 
the Iberian Peninsula and in Turkey. Here, two or three 
endemic related lineages, respectively, occupy geo- 
graphically different regions that likely acted as distinct 
refugia within a main refugium. It is noteworthy that in 
Iberia, where the species currently has a discontinuous 
distribution, genetic discontinuities do not fully corres- 
pond to geographic isolates, suggesting recent fragmen- 
tation dynamics. 

A scenario of survival in multiple refugia within a 
southern peninsula has been inferred for many species 
of the continental Europe fauna and detailed studies on 
amphibians and reptiles provide the best evidence (see 
[57] for a recent review). For example, there is evidence 
for multiple Iberian refugia in Salamandra salamandra 
[58], multiple Italian refugia in Triturus carnifex [59] 
and various Balkan refugia in Lacerta viridis [60]. To 
our knowledge, however, P. muralis is the first reported 
case of differentiation into multiple refugia within all 
Mediterranean peninsulas. In fact, within P. muralis, di- 
versity and divergence levels are not strikingly different 
among these regions (Table 2), suggesting that all Medi- 
terranean peninsulas significantly and simultaneously 
contributed to the long-term persistence of this species 
throughout the Pleistocene glaciations. Therefore, the 
case of P. muralis not only conforms to the scenario of 
"refugia within refugia" [56], but it also sets a new 
phylogeographic pattern of "refugia within all refugia". 

According to the high number of lineages found in 
southern European peninsulas, the phylogeographic pat- 
tern of P. muralis fits the classical southern refugia 
model to some extent. However, while for many species 
of continental Europe these Mediterranean peninsulas 
were both glacial refugia and source areas for northward 
postglacial colonization [2,4], in the case of P. muralis 
peninsular lineages are not found outside of these areas. 
This means that southern peninsulas acted as glacial re- 
fugia but not as sources for postglacial expansion. A 
possible explanation for such an unusual pattern is that 
northern areas outside Mediterranean peninsulas were 
already occupied by the species. In fact, various distinct 
lineages of P, muralis were found outside southern pen- 
insulas, in the Pyrenees and France (clades 1 and 6), 
Northern Italy (clade 5), Slovenia and Austria (clade 11) 
and the Central Balkans (clade 4), reaching the highest 
latitudes of the species range. Phylogenetic analyses and 
permution tests showed that these lineages are well 
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differentiated and evolutionarily independent from south- 
ern lineages. The occurrence of two allopatric lineages in 
Western and Eastern France (clades 1 and 6), which are 
unrelated and considerably divergent from each other and 
relative to lineages occurring in Iberia (clade 2 and 3), and 
on the other side of the Alpine arch in North Italy (clade 
5), provides evidence for a long-term persistence of P. 
muralis in these areas during Pleistocene glaciations, 
within separate refugia. Considering that the highest gen- 
etic diversity of these lineages is found along the southern 
portion of their range in SW and SE France respectively 
(data not shown), we can hypothesize that the refugia 
where these lineages differentiated were located in prox- 
imity of the Mediterranean coast and acted as source areas 
for post-glacial recolonization of Northern France and 
Germany. This hypothesis has been tested and validated 
by a recent study focused on P. muralis populations from 
France and Luxembourg [61], and it is further supported 
by palaeoclimatic and phylogeographic studies which indi- 
cate southern France as an area of prolonged ecological 
stability where climatic oscillations were attenuated [62] 
and where genetic signatures of refugia have been found 
among multiple species, including plants and animals (e.g. 
[63,64]). A similar reasoning applies to the allopatric mito- 
chondrial lineages found in Northern Italy, Northern 
Adriatic, Eastern Slovenia and Austria, and Central 
Balkans (clades 5, 21, 11, and 4), which are allopatric and 
reciprocally monophyletic, thus suggesting long-term 
isolation in multiple ice-age refugia. The long-term per- 
sistence of temperate species in Northern Italy (includ- 
ing the Po Plain and the southern foothills of the Alps) 
and in the Northern Adriatic has been inferred for many 
amphibians and reptiles in recent phylogeographic stud- 
ies based on the occurrence of divergent lineages en- 
demic to these areas (e.g., Hyla intermedia, Pelobates 
fuscus, Pelophylax lessonae, Triturus carnifex, Podarcis 
sicula; see [59,65]). Likewise, in the area spanning from 
the Eastern Alps to the South- Western Pannonian re- 
gion, palaeobotanical, paleoclimatic, and genetic data 
indicate that temperate species such as the tree Fagus 
sylvatica and the butterfly Erebia medusa survived in 
isolated refugia during glacial periods [20,21,66,67]. Fi- 
nally, a refugial area in the central Balkans and south- 
western Carpathians has already been suggested for 
many species showing high genetic diversity and distinct 
lineages in this area where also broad-leaved trees sur- 
vived glacial periods (e.g. [16,66-68]). 

Altogether these phylogeographic, palaeobotanical 
and palaeoclimatic studies provide evidence for the 
long-term persistence of temperate organisms outside 
the traditional refugia in southern peninsulas, and 
underline the prominent contribution of these northern 
refugia to the present-day genetic pools of many tem- 
perate species. 



Conclusions 

In recent decades, the phylogeography of the Western 
Palaearctic has been mainly centered around the south- 
ern refugia model, which sets a sharp dichotomy be- 
tween Mediterranean peninsulas and Northern Europe, 
the former being areas of long-term persistence of tem- 
perate species and sources for the postglacial assembly 
of northern biotas. Our data do not strictly conform to 
this model: a pattern of northern purity was not detected, 
whereas several evolutionary units endemic to northern 
regions were identifed. Overall, the occurrence of many 
reciprocally monophyletic lineages within P. muralis sug- 
gests a history of long-term persistence in multiple glacial 
refugia across its distribution, providing a paradigm of ice- 
age survival of temperate species in Mediterranean and 
extra-Mediterranean refugia [12]. Contrary to previous hy- 
potheses [25,27,28] this species did not occupy most of its 
present range only recently, either from the Italian Penin- 
sula or from any other southern peninsular refugium. 

In the last two decades mounting evidence for glacial 
survival in northern refugia through the Pleistocene has 
been found among various organisms representing dif- 
ferent biogeographical groups ([12] for a review). These 
studies suggest that the assembly of northern European 
biota was not simply the outcome of long-distance post- 
glacial re-colonization from southern peninsulas [6,7]. 
Rather, extra-Mediterranean refugia allowed the long- 
term persistence and differentiation of distinct genetic 
lineages of many organisms in the Western Palaearctic 
and acted as source areas for post-glacial colonization of 
Central and Northern Europe [12]. These biogeographic 
insights in the Western Palaearctic are paralleled by 
phylogeographic and palaeobotanical evidence for survival 
during ice ages of populations of temperate organisms in 
northern refugia located in the Eastern Palaearctic [69], as 
well as in the Nearctic region [70-72] providing directions 
for future research in the Holarctic ecozone concerning 
the role of northern regions as areas of long-term persist- 
ence of temperate biotas. 

Understanding the relative contribution of Mediterra- 
nean and extra-Mediterranean refugia to the current 
patterns of genetic diversity may allow us to significantly 
improve our current knowledge regarding glacial refugia 
dynamics, the extent and the routes of dispersal between 
subregions and their consequences in terms of evolu- 
tionary and demographic processes [9-12,73]. Moreover, 
our understanding of northern refugial biodiversity has 
deep consequences for conservation. The survival of tem- 
perate biota in northern refugia across a matrix of unsuit- 
able glacial landscapes challenges the view that organisms 
must have gone extinct in northern regions during glacia- 
tions, and suggests that, in many cases, postglacial colon- 
isation routes - and thus migration capacity - tracking the 
expansion of suitable habitats have been much smaller 
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than previously thought [9-12,68,74]. Our appraisal of past 
species resilience, and of their ability to track suitable hab- 
itats in response to past climatic shifts, is crucial to better 
forecast their future performance under global climate 
change scenarios [10,73,75-77], and thus essential for 
more informed decisions aimed at long-term conservation 
of biodiversity. 
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